Explore GLMM models for DNA methylation data
KE Lotterhos
9/26/2019
Install packages
packages_needed <- c("lme4", "multcomp", "devtools", "car", "MCMCglmm", "brms", "rstan", "emmeans", "tidybayes")
for (i in 1:length(packages_needed)){
if(!(packages_needed[i] %in% installed.packages())){install.packages(packages_needed[i])}
}
for (i in 1:length(packages_needed)){
library( packages_needed[i], character.only = TRUE)
}## Loading required package: Matrix## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Loading required package: MASS##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyser## Loading required package: usethis## Loading required package: carData## Registered S3 methods overwritten by 'car':
## method from
## influence.merMod lme4
## cooks.distance.influence.merMod lme4
## dfbeta.influence.merMod lme4
## dfbetas.influence.merMod lme4## Loading required package: coda## Loading required package: ape## Loading required package: Rcpp## Registered S3 method overwritten by 'xts':
## method from
## as.zoo.xts zoo## Loading 'brms' package (version 2.10.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following object is masked from 'package:MCMCglmm':
##
## me## The following object is masked from 'package:survival':
##
## kidney## The following object is masked from 'package:lme4':
##
## ngrps## Loading required package: StanHeaders## Loading required package: ggplot2## rstan (Version 2.19.2, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)##
## Attaching package: 'rstan'## The following object is masked from 'package:coda':
##
## traceplot##
## Attaching package: 'emmeans'## The following object is masked from 'package:devtools':
##
## testLoad data
Katie’s working directory for this markdown: setwd(“~/Documents/GitHub/2017OAExp_Oysters/markdown_files/DNAm/exploreBinomModels”)
You may have to set a different one.
## 17005 17007 17013 17019 17069 17070 17072 17079 17090 17094 17108 17122
## 1: 36 56 45 71 57 36 54 12 64 54 36 34
## 2: 32 57 43 75 47 36 51 0 42 45 0 35
## 3: 60 95 62 112 98 57 69 24 126 115 58 41
## 4: 39 57 49 59 59 33 38 19 71 65 36 24
## 5: 43 73 51 67 71 37 42 21 76 80 38 33
## 6: 26 15 23 9 18 6 19 27 17 13 13 16
## 17130 17142 17145 17162 17174 17176 17178 17181 17203 17211 17213
## 1: 28 47 55 51 58 76 77 26 65 61 35
## 2: 16 34 60 56 35 44 82 7 77 41 21
## 3: 56 64 97 81 77 107 84 68 105 94 65
## 4: 41 45 69 54 68 60 50 37 65 52 50
## 5: 39 48 64 56 61 79 58 47 74 59 45
## 6: 14 14 24 25 21 24 16 35 34 21 8## 17005 17007 17013 17019 17069 17070 17072 17079 17090 17094 17108 17122
## 1: 0 4 3 4 5 3 3 3 20 6 5 1
## 2: 8 3 11 5 20 9 11 16 46 26 46 2
## 3: 5 2 5 6 5 7 5 2 3 6 3 14
## 4: 7 13 7 10 14 10 3 2 9 23 5 11
## 5: 3 1 5 3 6 3 0 0 3 9 3 3
## 6: 14 6 0 4 0 0 0 1 6 1 1 3
## 17130 17142 17145 17162 17174 17176 17178 17181 17203 17211 17213
## 1: 2 1 3 2 4 0 3 18 2 4 6
## 2: 16 20 3 1 28 40 1 40 1 29 22
## 3: 2 3 7 4 18 7 10 2 5 2 2
## 4: 7 5 11 7 7 20 14 9 13 16 4
## 5: 5 4 10 3 12 3 8 1 4 3 7
## 6: 6 4 0 0 2 3 22 2 1 8 0## chr_1 pos cg_strand motif tri cg_index match_num gene_index
## 5449 NC_035780.1 100576 + CG CGG 5449 1 5
## 5451 NC_035780.1 100582 + CG CGA 5451 1 5
## 5453 NC_035780.1 100635 + CG CGC 5453 1 5
## 5455 NC_035780.1 100644 + CG CGG 5455 1 5
## 5457 NC_035780.1 100652 + CG CGG 5457 1 5
## 12986 NC_035780.1 258624 - CG CGA 12986 1 18
## chr_2 method feature start end score gene_strand phase
## 5449 NC_035780.1 Gnomon gene 99840 106460 . + .
## 5451 NC_035780.1 Gnomon gene 99840 106460 . + .
## 5453 NC_035780.1 Gnomon gene 99840 106460 . + .
## 5455 NC_035780.1 Gnomon gene 99840 106460 . + .
## 5457 NC_035780.1 Gnomon gene 99840 106460 . + .
## 12986 NC_035780.1 Gnomon gene 258108 272839 . - .
## gene_ID gene_id
## 5449 gene4 LOC111120752
## 5451 gene4 LOC111120752
## 5453 gene4 LOC111120752
## 5455 gene4 LOC111120752
## 5457 gene4 LOC111120752
## 12986 gene17 LOC111124802## sample_name sample_index treatment timepoint Day D9_C D9_E D80_C D80_E
## 1 RNA17005 1 400 6 80 0 0 1 0
## 2 RNA17007 2 400 6 80 0 0 1 0
## 3 RNA17013 3 400 3 9 1 0 0 0
## 4 RNA17019 4 2800 6 80 0 0 0 1
## 5 RNA17069 5 400 3 9 1 0 0 0
## 6 RNA17070 6 400 3 9 1 0 0 0
## shelf tank tankID population extraction_order seq_order lane
## 1 1 1 1 2 21 22 2
## 2 1 2 2 1 17 10 1
## 3 2 2 5 1 11 5 1
## 4 4 3 9 2 14 20 2
## 5 1 1 1 3 3 1 1
## 6 1 2 2 1 7 13 2
## gw_tapestation_RINe gw_tapestation_conc read_num dry_wgtcorr
## 1 8.3 50.2 40080928 0.715
## 2 9.0 75.9 42806943 0.177
## 3 8.7 46.1 42332976 1.110
## 4 8.9 61.4 45046863 0.707
## 5 8.4 41.3 42702450 1.040
## 6 9.5 71.9 41410387 1.250
## HN_Condition Simple_Condition epf_pH diff_pH env_pH Treatment
## 1 1.529600 23.983726 6.846758 -0.858465531 7.705224 400
## 2 7.687998 8.656369 7.455280 -0.381151365 7.836431 400
## 3 2.124024 6.307993 7.385796 -0.434471881 7.820268 400
## 4 2.709414 24.715950 7.149452 0.088113528 7.061338 2800
## 5 1.390289 7.889723 7.817657 0.008861369 7.808795 400
## 6 3.181037 4.404745 6.786949 -1.016110687 7.803059 400
## Time Pop Lane SFV ID
## 1 80 2 2 80.400 17005
## 2 80 1 1 80.400 17007
## 3 09 1 1 09.400 17013
## 4 80 2 2 80.2800 17019
## 5 09 3 1 09.400 17069
## 6 09 1 2 09.400 17070## [1] TRUE## [1] TRUECreate some interesting cases
# real data
a <- data.frame(mc=as.numeric(mC[1,]),
uc=as.numeric(uC[1,]),
level=meta_samp$level,
Treatment = meta_samp$Treatment,
Time = meta_samp$Time,
tankID = meta_samp$tankID,
Pop = meta_samp$Pop,
size = as.numeric(mC[1,])+as.numeric(uC[1,]))
a$level <- factor(a$level, levels=c("400_09", "400_80", "2800_09", "2800_80"))
levels(a$level)## [1] "400_09" "400_80" "2800_09" "2800_80"# real data with all 0 for unmeth in one treatment
b <- a
b$uc[b$level=='400_09'] <- 0
b$size <- b$mc+b$uc
b## mc uc level Treatment Time tankID Pop size
## 1 36 0 400_80 400 80 1 2 36
## 2 56 4 400_80 400 80 2 1 60
## 3 45 0 400_09 400 09 5 1 45
## 4 71 4 2800_80 2800 80 9 2 75
## 5 57 0 400_09 400 09 1 3 57
## 6 36 0 400_09 400 09 2 1 36
## 7 54 3 2800_09 2800 09 9 2 57
## 8 12 3 400_80 400 80 6 2 15
## 9 64 20 2800_09 2800 09 11 1 84
## 10 54 6 2800_80 2800 80 11 3 60
## 11 36 5 2800_09 2800 09 12 3 41
## 12 34 1 2800_09 2800 09 7 3 35
## 13 28 2 2800_80 2800 80 12 2 30
## 14 47 1 2800_09 2800 09 10 2 48
## 15 55 3 2800_80 2800 80 10 3 58
## 16 51 0 400_09 400 09 6 2 51
## 17 58 0 400_09 400 09 3 2 58
## 18 76 0 400_09 400 09 4 3 76
## 19 77 3 2800_80 2800 80 8 1 80
## 20 26 18 2800_09 2800 09 8 1 44
## 21 65 2 400_80 400 80 4 3 67
## 22 61 4 400_80 400 80 5 1 65
## 23 35 6 2800_80 2800 80 7 3 41# real data with all 0 for meth in one treatment
c <- a
c$mc[c$level=='400_09'] <- 0
c[18,2] <-5 # so no 0's in both columns
c$size <- c$mc+c$uc
# real data with all one count for meth in one treatment
# compare to c for zero variance
d <- a
d$mc[d$level=='400_09'] <- 1
d$uc[d$level=='400_09'] <- 20
d$size <- d$mc+d$uc
d## mc uc level Treatment Time tankID Pop size
## 1 36 0 400_80 400 80 1 2 36
## 2 56 4 400_80 400 80 2 1 60
## 3 1 20 400_09 400 09 5 1 21
## 4 71 4 2800_80 2800 80 9 2 75
## 5 1 20 400_09 400 09 1 3 21
## 6 1 20 400_09 400 09 2 1 21
## 7 54 3 2800_09 2800 09 9 2 57
## 8 12 3 400_80 400 80 6 2 15
## 9 64 20 2800_09 2800 09 11 1 84
## 10 54 6 2800_80 2800 80 11 3 60
## 11 36 5 2800_09 2800 09 12 3 41
## 12 34 1 2800_09 2800 09 7 3 35
## 13 28 2 2800_80 2800 80 12 2 30
## 14 47 1 2800_09 2800 09 10 2 48
## 15 55 3 2800_80 2800 80 10 3 58
## 16 1 20 400_09 400 09 6 2 21
## 17 1 20 400_09 400 09 3 2 21
## 18 1 20 400_09 400 09 4 3 21
## 19 77 3 2800_80 2800 80 8 1 80
## 20 26 18 2800_09 2800 09 8 1 44
## 21 65 2 400_80 400 80 4 3 67
## 22 61 4 400_80 400 80 5 1 65
## 23 35 6 2800_80 2800 80 7 3 41Test different models for real data
boxplot(mc/(mc+uc) ~ Treatment*Time, data=a, las=2)
# No random effect of tank
m1 <- glm(cbind(mc,uc)~ Treatment*Time, data=a,family = "binomial")
summary(m1)##
## Call:
## glm(formula = cbind(mc, uc) ~ Treatment * Time, family = "binomial",
## data = a)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.0350 -0.7922 0.0671 1.0852 3.1420
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.94444 0.24884 11.833 < 2e-16 ***
## Treatment2800 -1.25112 0.29425 -4.252 2.12e-05 ***
## Time80 -0.07131 0.37841 -0.188 0.8505
## Treatment2800:Time80 0.96826 0.46114 2.100 0.0358 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 94.020 on 22 degrees of freedom
## Residual deviance: 67.161 on 19 degrees of freedom
## AIC: 139.65
##
## Number of Fisher Scoring iterations: 5boxplot(m1$residuals~a$tankID)
boxplot(m1$residuals~a$Pop)
# definitely some tank effects in the residuals
vif(m1) # variance inflation factor## Treatment Time Treatment:Time
## 1.696074 3.078354 3.586870AIC(m1)## [1] 139.6509m1.0 <- glm(cbind(mc,uc)~ Treatment*Time, data=a,family = "quasibinomial")
summary(m1.0)##
## Call:
## glm(formula = cbind(mc, uc) ~ Treatment * Time, family = "quasibinomial",
## data = a)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.0350 -0.7922 0.0671 1.0852 3.1420
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.94444 0.45704 6.442 3.55e-06 ***
## Treatment2800 -1.25112 0.54045 -2.315 0.032 *
## Time80 -0.07131 0.69502 -0.103 0.919
## Treatment2800:Time80 0.96826 0.84697 1.143 0.267
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 3.373489)
##
## Null deviance: 94.020 on 22 degrees of freedom
## Residual deviance: 67.161 on 19 degrees of freedom
## AIC: NA
##
## Number of Fisher Scoring iterations: 5boxplot(m1.0$residuals~a$tankID)
boxplot(m1.0$residuals~a$Pop)
# definitely some tank effects in the residuals
vif(m1.0) # variance inflation factor## Treatment Time Treatment:Time
## 1.696074 3.078354 3.586870m1.1 <- glmer(cbind(mc,uc)~ Treatment*Time + (1|tankID), data=a,family = "binomial")
summary(m1.1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(mc, uc) ~ Treatment * Time + (1 | tankID)
## Data: a
##
## AIC BIC logLik deviance df.resid
## 132.9 138.6 -61.5 122.9 18
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9594 -0.5257 -0.1039 0.9119 1.9449
##
## Random effects:
## Groups Name Variance Std.Dev.
## tankID (Intercept) 0.1991 0.4462
## Number of obs: 23, groups: tankID, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.97896 0.30990 9.612 < 2e-16 ***
## Treatment2800 -1.16124 0.39613 -2.931 0.00337 **
## Time80 -0.09855 0.39200 -0.251 0.80150
## Treatment2800:Time80 1.00737 0.47517 2.120 0.03400 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tr2800 Time80
## Tretmnt2800 -0.775
## Time80 -0.524 0.407
## Trt2800:T80 0.433 -0.471 -0.825boxplot(residuals(m1.1)~a$tankID)
boxplot(residuals(m1.1)~a$Pop)
# this seems to be better, but still not perfectly accounting for residuals across tanks
vif(m1.1) # the inflation in size of the confidence ellipse or ellipsoid for the coefficients of the term in comparison with what would be obtained for orthogonal data.## Treatment Time Treatment:Time
## 1.286612 3.142746 3.370736AIC(m1.1)## [1] 132.9438 # improvement in model fit compared to m0Test different models for fake data “b”
One treatment has 100% methylation in all individuals
boxplot(mc/(mc+uc) ~ Treatment*Time, data=b, las=2)
# No random effect of tank
m1 <- glm(cbind(mc,uc)~ Treatment*Time, data=b,family = "binomial")
summary(m1)##
## Call:
## glm(formula = cbind(mc, uc) ~ Treatment * Time, family = "binomial",
## data = b)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.0350 -0.3595 0.0002 0.7722 3.1420
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 21.68 1720.93 0.013 0.990
## Treatment2800 -19.99 1720.93 -0.012 0.991
## Time80 -18.81 1720.93 -0.011 0.991
## Treatment2800:Time80 19.70 1720.93 0.011 0.991
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 131.486 on 22 degrees of freedom
## Residual deviance: 57.225 on 19 degrees of freedom
## AIC: 114.7
##
## Number of Fisher Scoring iterations: 17boxplot(m1$residuals~b$tankID)
boxplot(m1$residuals~b$Pop)
# definitely some tank effects in the residuals
vif(m1) # variance inflation factor## Treatment Time Treatment:Time
## 30476288 55314099 46482754 # hmmm
m1.1 <- glmer(cbind(mc,uc)~ Treatment*Time + (1|tankID), data=b,family = "binomial")## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.012937
## (tol = 0.001, component 1)summary(m1.1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(mc, uc) ~ Treatment * Time + (1 | tankID)
## Data: b
##
## AIC BIC logLik deviance df.resid
## 107.5 113.2 -48.8 97.5 18
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.82988 -0.39891 0.00005 0.34430 2.00964
##
## Random effects:
## Groups Name Variance Std.Dev.
## tankID (Intercept) 0.267 0.5167
## Number of obs: 23, groups: tankID, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 23.711322 0.005465 4339 <2e-16 ***
## Treatment2800 -21.876194 0.005459 -4007 <2e-16 ***
## Time80 -20.794207 0.005465 -3805 <2e-16 ***
## Treatment2800:Time80 21.706373 0.005459 3976 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tr2800 Time80
## Tretmnt2800 0.000
## Time80 0.001 0.000
## Trt2800:T80 0.000 0.000 0.000
## convergence code: 0
## Model failed to converge with max|grad| = 0.012937 (tol = 0.001, component 1)boxplot(residuals(m1.1)~b$tankID)
boxplot(residuals(m1.1)~b$Pop)
# this seems to be better, but we get that strange error message
vif(m1.1) # the inflation in size of the confidence ellipse or ellipsoid for the coefficients of the term in comparison with what would be obtained for orthogonal data.## Treatment Time Treatment:Time
## 1 1 1Error: “Model failed to converge with max|grad| = 0.012937 (tol = 0.001, component 1)”
Test different models for fake data “c”
One treatment has 0% methylation in all individuals
boxplot(mc/(mc+uc) ~ Treatment*Time, data=c, las=2)
# No random effect of tank
m1 <- glm(cbind(mc,uc)~ Treatment*Time, data=c,family = "binomial")
summary(m1)##
## Call:
## glm(formula = cbind(mc, uc) ~ Treatment * Time, family = "binomial",
## data = c)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.0350 -0.3595 -0.0002 0.7722 3.1420
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -19.30 2006.20 -0.010 0.992
## Treatment2800 20.99 2006.20 0.010 0.992
## Time80 22.17 2006.20 0.011 0.991
## Treatment2800:Time80 -21.28 2006.20 -0.011 0.992
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 175.801 on 22 degrees of freedom
## Residual deviance: 57.225 on 19 degrees of freedom
## AIC: 114.7
##
## Number of Fisher Scoring iterations: 17boxplot(m1$residuals~c$tankID)
boxplot(m1$residuals~c$Pop)
# definitely some tank effects in the residuals
vif(m1) # variance inflation factor## Treatment Time Treatment:Time
## 41417750 75172723 63170786 # hmmm not good
m1.1 <- glmer(cbind(mc,uc)~ Treatment*Time + (1|tankID), data=c,family = "binomial")## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
## control$checkConv, : Model failed to converge with max|grad| = 0.0129364
## (tol = 0.001, component 1)summary(m1.1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(mc, uc) ~ Treatment * Time + (1 | tankID)
## Data: c
##
## AIC BIC logLik deviance df.resid
## 107.5 113.2 -48.8 97.5 18
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.82997 -0.39894 -0.00011 0.34429 2.00970
##
## Random effects:
## Groups Name Variance Std.Dev.
## tankID (Intercept) 0.267 0.5167
## Number of obs: 23, groups: tankID, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -20.057445 0.005466 -3670 <2e-16 ***
## Treatment2800 21.892588 0.005460 4010 <2e-16 ***
## Time80 22.974561 0.005465 4204 <2e-16 ***
## Treatment2800:Time80 -22.062452 0.005460 -4041 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tr2800 Time80
## Tretmnt2800 0.000
## Time80 0.001 0.000
## Trt2800:T80 0.000 0.000 0.000
## convergence code: 0
## Model failed to converge with max|grad| = 0.0129364 (tol = 0.001, component 1)boxplot(residuals(m1.1)~c$tankID)
boxplot(residuals(m1.1)~c$Pop)
# this seems to be better, but we get that strange error message
vif(m1.1) # the inflation in size of the confidence ellipse or ellipsoid for the coefficients of the term in comparison with what would be obtained for orthogonal data.## Treatment Time Treatment:Time
## 1 1 1Error: “Model failed to converge with max|grad| = 0.0129364 (tol = 0.001, component 1)”
Test different models for fake data “d”
One treatment has a count of 1 methylation in all individuals. This can be compared to dataset “c” to show the problem comes from fixation of counts of 0 and not from zero variance in a treatment.
boxplot(mc/(mc+uc) ~ Treatment*Time, data=d, las=2)
# No random effect of tank
m1 <- glm(cbind(mc,uc)~ Treatment*Time, data=d,family = "binomial")
summary(m1)##
## Call:
## glm(formula = cbind(mc, uc) ~ Treatment * Time, family = "binomial",
## data = d)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.0350 -0.3595 0.0000 0.7722 3.1420
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.9957 0.4183 -7.161 8e-13 ***
## Treatment2800 4.6891 0.4468 10.494 <2e-16 ***
## Time80 5.8689 0.5062 11.593 <2e-16 ***
## Treatment2800:Time80 -4.9719 0.5707 -8.712 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 491.062 on 22 degrees of freedom
## Residual deviance: 57.225 on 19 degrees of freedom
## AIC: 126.41
##
## Number of Fisher Scoring iterations: 5boxplot(m1$residuals~d$tankID)
boxplot(m1$residuals~d$Pop)
# definitely some tank effects in the residuals
vif(m1) # variance inflation factor## Treatment Time Treatment:Time
## 2.796285 5.075224 5.264923 # interesting, much better behaved
m1.1 <- glmer(cbind(mc,uc)~ Treatment*Time + (1|tankID), data=d,family = "binomial")
# note no error message
summary(m1.1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(mc, uc) ~ Treatment * Time + (1 | tankID)
## Data: d
##
## AIC BIC logLik deviance df.resid
## 120.3 126.0 -55.2 110.3 18
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9463 -0.3925 -0.0720 0.4843 1.9514
##
## Random effects:
## Groups Name Variance Std.Dev.
## tankID (Intercept) 0.2047 0.4524
## Number of obs: 23, groups: tankID, 12
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.0704 0.4635 -6.624 3.51e-11 ***
## Treatment2800 4.8897 0.5314 9.201 < 2e-16 ***
## Time80 5.9682 0.5254 11.359 < 2e-16 ***
## Treatment2800:Time80 -5.0591 0.5892 -8.587 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Tr2800 Time80
## Tretmnt2800 -0.880
## Time80 -0.747 0.660
## Trt2800:T80 0.665 -0.669 -0.890boxplot(residuals(m1.1)~d$tankID)
boxplot(residuals(m1.1)~d$Pop)
# this seems to be better, but we get that strange error message
vif(m1.1) # the inflation in size of the confidence ellipse or ellipsoid for the coefficients of the term in comparison with what would be obtained for orthogonal data.## Treatment Time Treatment:Time
## 1.877977 5.008708 5.107818Better understand what is causing error in datasets b and c
These show the error is caused by all individuals within a treatment having 0% or 100% methylation. Clearly, these are cases that are very interesting!
Googling “Warning message: In checkConv(attr(opt,”derivs“), opt\(par, ctrl = control\)checkConv, : Model failed to converge with max|grad| = 0.0129364 (tol = 0.001, component 1)”
This seems to be a very common error message, but I couldn’t find a specific post that highlights the problem we discovered here.
lme4 is based on restricted maximum likelihood, which appears to have problems converging when all individuals within a treatment have 0% or 100% methylation
Note: Alan also tried different converging algorithms and also had problems
Bayesian models offer a good alternative…
Let’s try Bayesian MCMCglmm!
- Read Jared H’s book in the Lotterhos Lab Reading / stats
- Good example here
# These are default priors I found in some examples. I need to read up on these and
prior1=list(R=list(V=1, nu=0.002), G=list(G1=list(V=1, nu=0.002)))
prior2=list(R=list(V=1, nu=0.002), G=list(G1=list(V=1, nu=1,
alpha.mu=0, alpha.V=100)))
# Real data
m2 <- MCMCglmm(cbind(mc,uc)~ Treatment*Time, random= ~tankID,
data=a,family = "multinomial2", prior=prior1)##
## MCMC iteration = 0
##
## Acceptance ratio for liability set 1 = 0.000522
##
## MCMC iteration = 1000
##
## Acceptance ratio for liability set 1 = 0.443000
##
## MCMC iteration = 2000
##
## Acceptance ratio for liability set 1 = 0.443826
##
## MCMC iteration = 3000
##
## Acceptance ratio for liability set 1 = 0.442348
##
## MCMC iteration = 4000
##
## Acceptance ratio for liability set 1 = 0.422348
##
## MCMC iteration = 5000
##
## Acceptance ratio for liability set 1 = 0.421435
##
## MCMC iteration = 6000
##
## Acceptance ratio for liability set 1 = 0.415739
##
## MCMC iteration = 7000
##
## Acceptance ratio for liability set 1 = 0.428391
##
## MCMC iteration = 8000
##
## Acceptance ratio for liability set 1 = 0.434130
##
## MCMC iteration = 9000
##
## Acceptance ratio for liability set 1 = 0.433870
##
## MCMC iteration = 10000
##
## Acceptance ratio for liability set 1 = 0.417217
##
## MCMC iteration = 11000
##
## Acceptance ratio for liability set 1 = 0.426478
##
## MCMC iteration = 12000
##
## Acceptance ratio for liability set 1 = 0.426739
##
## MCMC iteration = 13000
##
## Acceptance ratio for liability set 1 = 0.424826plot(cbind(m2$VCV), pch=19, cex=0.2)
par(mfrow=c(4,2), mar=c(2,2,1,0))
plot(m2$Sol, auto.layout=F)
# no autocorrelation in the traces, good!
summary(m2)##
## Iterations = 3001:12991
## Thinning interval = 10
## Sample size = 1000
##
## DIC: 652.7295
##
## G-structure: ~tankID
##
## post.mean l-95% CI u-95% CI eff.samp
## tankID 0.1182 0.0004026 0.5282 333.7
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.734 0.1008 1.596 641.3
##
## Location effects: cbind(mc, uc) ~ Treatment * Time
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 3.10206 2.20106 3.98942 1000 <0.001 ***
## Treatment2800 -1.00286 -2.37882 0.05296 1000 0.098 .
## Time80 -0.11617 -1.29458 1.35864 1126 0.852
## Treatment2800:Time80 0.69285 -0.99875 2.39050 1000 0.400
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # This result is interesting, more conservative P-values than
# what we got with the glmer
# Fake data with 0's for methylation
m2 <- MCMCglmm(cbind(mc,uc)~ Treatment*Time, random= ~tankID,
data=c,family = "multinomial2", prior=prior1)##
## MCMC iteration = 0
##
## Acceptance ratio for liability set 1 = 0.000696
##
## MCMC iteration = 1000
##
## Acceptance ratio for liability set 1 = 0.406609
##
## MCMC iteration = 2000
##
## Acceptance ratio for liability set 1 = 0.411304
##
## MCMC iteration = 3000
##
## Acceptance ratio for liability set 1 = 0.425304
##
## MCMC iteration = 4000
##
## Acceptance ratio for liability set 1 = 0.385130
##
## MCMC iteration = 5000
##
## Acceptance ratio for liability set 1 = 0.379652
##
## MCMC iteration = 6000
##
## Acceptance ratio for liability set 1 = 0.378000
##
## MCMC iteration = 7000
##
## Acceptance ratio for liability set 1 = 0.363391
##
## MCMC iteration = 8000
##
## Acceptance ratio for liability set 1 = 0.373913
##
## MCMC iteration = 9000
##
## Acceptance ratio for liability set 1 = 0.364304
##
## MCMC iteration = 10000
##
## Acceptance ratio for liability set 1 = 0.381913
##
## MCMC iteration = 11000
##
## Acceptance ratio for liability set 1 = 0.397522
##
## MCMC iteration = 12000
##
## Acceptance ratio for liability set 1 = 0.362913
##
## MCMC iteration = 13000
##
## Acceptance ratio for liability set 1 = 0.382174plot(cbind(m2$VCV), pch=19, cex=0.2)
par(mfrow=c(4,2), mar=c(2,2,1,0))
plot(m2$Sol, auto.layout=F)
# some autocorrelation in the traces, not good!
summary(m2)##
## Iterations = 3001:12991
## Thinning interval = 10
## Sample size = 1000
##
## DIC: 516.6384
##
## G-structure: ~tankID
##
## post.mean l-95% CI u-95% CI eff.samp
## tankID 0.1811 0.0002463 0.8145 329.6
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.8991 0.1175 2.101 588.8
##
## Location effects: cbind(mc, uc) ~ Treatment * Time
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -17.467 -36.471 -3.973 3.675 <0.001 ***
## Treatment2800 19.611 5.817 38.573 4.158 <0.001 ***
## Time80 20.494 7.468 40.336 3.880 <0.001 ***
## Treatment2800:Time80 -19.949 -39.395 -6.386 4.094 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Fake data with 0's for unmenthlated
m2 <- MCMCglmm(cbind(mc,uc)~ Treatment*Time, random= ~tankID,
data=b,family = "multinomial2", prior=prior1)##
## MCMC iteration = 0
##
## Acceptance ratio for liability set 1 = 0.000478
##
## MCMC iteration = 1000
##
## Acceptance ratio for liability set 1 = 0.417261
##
## MCMC iteration = 2000
##
## Acceptance ratio for liability set 1 = 0.436217
##
## MCMC iteration = 3000
##
## Acceptance ratio for liability set 1 = 0.446130
##
## MCMC iteration = 4000
##
## Acceptance ratio for liability set 1 = 0.433000
##
## MCMC iteration = 5000
##
## Acceptance ratio for liability set 1 = 0.433609
##
## MCMC iteration = 6000
##
## Acceptance ratio for liability set 1 = 0.427304
##
## MCMC iteration = 7000
##
## Acceptance ratio for liability set 1 = 0.433783
##
## MCMC iteration = 8000
##
## Acceptance ratio for liability set 1 = 0.438174
##
## MCMC iteration = 9000
##
## Acceptance ratio for liability set 1 = 0.443826
##
## MCMC iteration = 10000
##
## Acceptance ratio for liability set 1 = 0.452696
##
## MCMC iteration = 11000
##
## Acceptance ratio for liability set 1 = 0.428435
##
## MCMC iteration = 12000
##
## Acceptance ratio for liability set 1 = 0.441609
##
## MCMC iteration = 13000
##
## Acceptance ratio for liability set 1 = 0.435913plot(cbind(m2$VCV), pch=19, cex=0.2)
par(mfrow=c(4,2), mar=c(2,2,1,0))
plot(m2$Sol, auto.layout=F)
# lots of autocorrelation in the traces, not good! Will need to look into convergence issues
summary(m2)##
## Iterations = 3001:12991
## Thinning interval = 10
## Sample size = 1000
##
## DIC: 516.414
##
## G-structure: ~tankID
##
## post.mean l-95% CI u-95% CI eff.samp
## tankID 0.1266 0.0002649 0.6104 386
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.9213 0.09074 2.021 526.7
##
## Location effects: cbind(mc, uc) ~ Treatment * Time
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) 22.156 10.388 33.429 5.912 <0.001 ***
## Treatment2800 -20.030 -31.375 -8.207 6.224 <0.001 ***
## Time80 -19.125 -30.787 -7.490 6.455 <0.001 ***
## Treatment2800:Time80 19.687 8.219 31.655 6.921 <0.001 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# let's see if we can do better by changing the prior and increasing
# the number of iterations
m2 <- MCMCglmm(cbind(mc,uc)~ Treatment*Time, random= ~tankID,
data=c,family = "multinomial2", prior=prior2,
thin = 20,
burnin = 3000,
nitt = 100000)##
## MCMC iteration = 0
##
## Acceptance ratio for liability set 1 = 0.000652
##
## MCMC iteration = 1000
##
## Acceptance ratio for liability set 1 = 0.407870
##
## MCMC iteration = 2000
##
## Acceptance ratio for liability set 1 = 0.418696
##
## MCMC iteration = 3000
##
## Acceptance ratio for liability set 1 = 0.387174
##
## MCMC iteration = 4000
##
## Acceptance ratio for liability set 1 = 0.377739
##
## MCMC iteration = 5000
##
## Acceptance ratio for liability set 1 = 0.379391
##
## MCMC iteration = 6000
##
## Acceptance ratio for liability set 1 = 0.354348
##
## MCMC iteration = 7000
##
## Acceptance ratio for liability set 1 = 0.386130
##
## MCMC iteration = 8000
##
## Acceptance ratio for liability set 1 = 0.408043
##
## MCMC iteration = 9000
##
## Acceptance ratio for liability set 1 = 0.374870
##
## MCMC iteration = 10000
##
## Acceptance ratio for liability set 1 = 0.384000
##
## MCMC iteration = 11000
##
## Acceptance ratio for liability set 1 = 0.370522
##
## MCMC iteration = 12000
##
## Acceptance ratio for liability set 1 = 0.369826
##
## MCMC iteration = 13000
##
## Acceptance ratio for liability set 1 = 0.379000
##
## MCMC iteration = 14000
##
## Acceptance ratio for liability set 1 = 0.385609
##
## MCMC iteration = 15000
##
## Acceptance ratio for liability set 1 = 0.346783
##
## MCMC iteration = 16000
##
## Acceptance ratio for liability set 1 = 0.374826
##
## MCMC iteration = 17000
##
## Acceptance ratio for liability set 1 = 0.387826
##
## MCMC iteration = 18000
##
## Acceptance ratio for liability set 1 = 0.387565
##
## MCMC iteration = 19000
##
## Acceptance ratio for liability set 1 = 0.388478
##
## MCMC iteration = 20000
##
## Acceptance ratio for liability set 1 = 0.391870
##
## MCMC iteration = 21000
##
## Acceptance ratio for liability set 1 = 0.379609
##
## MCMC iteration = 22000
##
## Acceptance ratio for liability set 1 = 0.383087
##
## MCMC iteration = 23000
##
## Acceptance ratio for liability set 1 = 0.363609
##
## MCMC iteration = 24000
##
## Acceptance ratio for liability set 1 = 0.375522
##
## MCMC iteration = 25000
##
## Acceptance ratio for liability set 1 = 0.376304
##
## MCMC iteration = 26000
##
## Acceptance ratio for liability set 1 = 0.373696
##
## MCMC iteration = 27000
##
## Acceptance ratio for liability set 1 = 0.384174
##
## MCMC iteration = 28000
##
## Acceptance ratio for liability set 1 = 0.377783
##
## MCMC iteration = 29000
##
## Acceptance ratio for liability set 1 = 0.379087
##
## MCMC iteration = 30000
##
## Acceptance ratio for liability set 1 = 0.385826
##
## MCMC iteration = 31000
##
## Acceptance ratio for liability set 1 = 0.385652
##
## MCMC iteration = 32000
##
## Acceptance ratio for liability set 1 = 0.370087
##
## MCMC iteration = 33000
##
## Acceptance ratio for liability set 1 = 0.381783
##
## MCMC iteration = 34000
##
## Acceptance ratio for liability set 1 = 0.384261
##
## MCMC iteration = 35000
##
## Acceptance ratio for liability set 1 = 0.365739
##
## MCMC iteration = 36000
##
## Acceptance ratio for liability set 1 = 0.393217
##
## MCMC iteration = 37000
##
## Acceptance ratio for liability set 1 = 0.394478
##
## MCMC iteration = 38000
##
## Acceptance ratio for liability set 1 = 0.385087
##
## MCMC iteration = 39000
##
## Acceptance ratio for liability set 1 = 0.384174
##
## MCMC iteration = 40000
##
## Acceptance ratio for liability set 1 = 0.368348
##
## MCMC iteration = 41000
##
## Acceptance ratio for liability set 1 = 0.374087
##
## MCMC iteration = 42000
##
## Acceptance ratio for liability set 1 = 0.379913
##
## MCMC iteration = 43000
##
## Acceptance ratio for liability set 1 = 0.370522
##
## MCMC iteration = 44000
##
## Acceptance ratio for liability set 1 = 0.388261
##
## MCMC iteration = 45000
##
## Acceptance ratio for liability set 1 = 0.383696
##
## MCMC iteration = 46000
##
## Acceptance ratio for liability set 1 = 0.379609
##
## MCMC iteration = 47000
##
## Acceptance ratio for liability set 1 = 0.373261
##
## MCMC iteration = 48000
##
## Acceptance ratio for liability set 1 = 0.371000
##
## MCMC iteration = 49000
##
## Acceptance ratio for liability set 1 = 0.387913
##
## MCMC iteration = 50000
##
## Acceptance ratio for liability set 1 = 0.370522
##
## MCMC iteration = 51000
##
## Acceptance ratio for liability set 1 = 0.367522
##
## MCMC iteration = 52000
##
## Acceptance ratio for liability set 1 = 0.381217
##
## MCMC iteration = 53000
##
## Acceptance ratio for liability set 1 = 0.376783
##
## MCMC iteration = 54000
##
## Acceptance ratio for liability set 1 = 0.369957
##
## MCMC iteration = 55000
##
## Acceptance ratio for liability set 1 = 0.382391
##
## MCMC iteration = 56000
##
## Acceptance ratio for liability set 1 = 0.391217
##
## MCMC iteration = 57000
##
## Acceptance ratio for liability set 1 = 0.398696
##
## MCMC iteration = 58000
##
## Acceptance ratio for liability set 1 = 0.366696
##
## MCMC iteration = 59000
##
## Acceptance ratio for liability set 1 = 0.366435
##
## MCMC iteration = 60000
##
## Acceptance ratio for liability set 1 = 0.379391
##
## MCMC iteration = 61000
##
## Acceptance ratio for liability set 1 = 0.385783
##
## MCMC iteration = 62000
##
## Acceptance ratio for liability set 1 = 0.375217
##
## MCMC iteration = 63000
##
## Acceptance ratio for liability set 1 = 0.364522
##
## MCMC iteration = 64000
##
## Acceptance ratio for liability set 1 = 0.345478
##
## MCMC iteration = 65000
##
## Acceptance ratio for liability set 1 = 0.381217
##
## MCMC iteration = 66000
##
## Acceptance ratio for liability set 1 = 0.397565
##
## MCMC iteration = 67000
##
## Acceptance ratio for liability set 1 = 0.374522
##
## MCMC iteration = 68000
##
## Acceptance ratio for liability set 1 = 0.381174
##
## MCMC iteration = 69000
##
## Acceptance ratio for liability set 1 = 0.368739
##
## MCMC iteration = 70000
##
## Acceptance ratio for liability set 1 = 0.370000
##
## MCMC iteration = 71000
##
## Acceptance ratio for liability set 1 = 0.379957
##
## MCMC iteration = 72000
##
## Acceptance ratio for liability set 1 = 0.371609
##
## MCMC iteration = 73000
##
## Acceptance ratio for liability set 1 = 0.383478
##
## MCMC iteration = 74000
##
## Acceptance ratio for liability set 1 = 0.383261
##
## MCMC iteration = 75000
##
## Acceptance ratio for liability set 1 = 0.371217
##
## MCMC iteration = 76000
##
## Acceptance ratio for liability set 1 = 0.376826
##
## MCMC iteration = 77000
##
## Acceptance ratio for liability set 1 = 0.376957
##
## MCMC iteration = 78000
##
## Acceptance ratio for liability set 1 = 0.378000
##
## MCMC iteration = 79000
##
## Acceptance ratio for liability set 1 = 0.349348
##
## MCMC iteration = 80000
##
## Acceptance ratio for liability set 1 = 0.383217
##
## MCMC iteration = 81000
##
## Acceptance ratio for liability set 1 = 0.395609
##
## MCMC iteration = 82000
##
## Acceptance ratio for liability set 1 = 0.382217
##
## MCMC iteration = 83000
##
## Acceptance ratio for liability set 1 = 0.372348
##
## MCMC iteration = 84000
##
## Acceptance ratio for liability set 1 = 0.372609
##
## MCMC iteration = 85000
##
## Acceptance ratio for liability set 1 = 0.372174
##
## MCMC iteration = 86000
##
## Acceptance ratio for liability set 1 = 0.382174
##
## MCMC iteration = 87000
##
## Acceptance ratio for liability set 1 = 0.381870
##
## MCMC iteration = 88000
##
## Acceptance ratio for liability set 1 = 0.377696
##
## MCMC iteration = 89000
##
## Acceptance ratio for liability set 1 = 0.380217
##
## MCMC iteration = 90000
##
## Acceptance ratio for liability set 1 = 0.377696
##
## MCMC iteration = 91000
##
## Acceptance ratio for liability set 1 = 0.379043
##
## MCMC iteration = 92000
##
## Acceptance ratio for liability set 1 = 0.363696
##
## MCMC iteration = 93000
##
## Acceptance ratio for liability set 1 = 0.378130
##
## MCMC iteration = 94000
##
## Acceptance ratio for liability set 1 = 0.377913
##
## MCMC iteration = 95000
##
## Acceptance ratio for liability set 1 = 0.388957
##
## MCMC iteration = 96000
##
## Acceptance ratio for liability set 1 = 0.390739
##
## MCMC iteration = 97000
##
## Acceptance ratio for liability set 1 = 0.368000
##
## MCMC iteration = 98000
##
## Acceptance ratio for liability set 1 = 0.374565
##
## MCMC iteration = 99000
##
## Acceptance ratio for liability set 1 = 0.372696
##
## MCMC iteration = 100000
##
## Acceptance ratio for liability set 1 = 0.365043 # further increasing the thinning results in worse density plots;
#
plot(cbind(m2$VCV), pch=19, cex=0.2)
par(mfrow=c(4,2), mar=c(2,2,1,0))
plot(m2$Sol, auto.layout=F)
summary(m2)##
## Iterations = 3001:99981
## Thinning interval = 20
## Sample size = 4850
##
## DIC: 517.0249
##
## G-structure: ~tankID
##
## post.mean l-95% CI u-95% CI eff.samp
## tankID 0.4497 1.271e-08 1.752 3477
##
## R-structure: ~units
##
## post.mean l-95% CI u-95% CI eff.samp
## units 0.8937 0.003763 2.094 3795
##
## Location effects: cbind(mc, uc) ~ Treatment * Time
##
## post.mean l-95% CI u-95% CI eff.samp pMCMC
## (Intercept) -122.85 -209.29 -18.81 1.483 <2e-04 ***
## Treatment2800 124.98 22.03 212.30 1.492 <2e-04 ***
## Time80 125.90 23.28 213.67 1.498 <2e-04 ***
## Treatment2800:Time80 -125.35 -212.30 -21.59 1.542 <2e-04 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # still some autocorrelation, but better than beforeA lot of autocorrelation in the trace implies that this is not converging!
In some cases we get an error message “In summary.glm(glm(cbind(MCMC_y, MCMC_y.additional) ~ 1, family =”quasibinomial“, : observations with zero weight not used for calculating dispersion”
This error also illustrates issues with the MCMC approach, as it relies on glms which are unhappy when all individuals within a treatment combination have 0% or 100% methylation.
Additionally, Burkner suggests that Metropolis-Hastings and Gibbs sampling is slow convergence for high-dimensional models with correlated parameters.
Let’s try Bayesian brms!
The github page has lots of great documentation
#vignette(package = "brms")
#help("brm")
# No random effect of tank on real data
m3 <- brm(mc | trials(size)~ Treatment*Time, data=a,family = binomial())## Compiling the C++ model## Start sampling##
## SAMPLING FOR MODEL '2f81065a9d1926dcf096558e0a7da943' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 2.9e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.29 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.068319 seconds (Warm-up)
## Chain 1: 0.067839 seconds (Sampling)
## Chain 1: 0.136158 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '2f81065a9d1926dcf096558e0a7da943' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.063504 seconds (Warm-up)
## Chain 2: 0.064401 seconds (Sampling)
## Chain 2: 0.127905 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '2f81065a9d1926dcf096558e0a7da943' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.7e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.06488 seconds (Warm-up)
## Chain 3: 0.059097 seconds (Sampling)
## Chain 3: 0.123977 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '2f81065a9d1926dcf096558e0a7da943' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.2e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.066593 seconds (Warm-up)
## Chain 4: 0.060986 seconds (Sampling)
## Chain 4: 0.127579 seconds (Total)
## Chain 4: # takes a few seconds
summary(m3)## Family: binomial
## Links: mu = logit
## Formula: mc | trials(size) ~ Treatment * Time
## Data: a (Number of observations: 23)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 2.97 0.25 2.51 3.49 1.00 1768
## Treatment2800 -1.27 0.29 -1.88 -0.72 1.00 1754
## Time80 -0.07 0.38 -0.81 0.69 1.00 1778
## Treatment2800:Time80 0.97 0.46 0.09 1.89 1.00 1658
## Tail_ESS
## Intercept 2028
## Treatment2800 1919
## Time80 2260
## Treatment2800:Time80 2251
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m3)
# gives a similar result ot the GLM
# Add the random effect of tank
m3.1 <- brm(mc | trials(size)~ Treatment*Time + (1|tankID), data=a,family = binomial())## Compiling the C++ model
## Start sampling##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 4e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.4 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.18596 seconds (Warm-up)
## Chain 1: 0.18513 seconds (Sampling)
## Chain 1: 0.37109 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 2.2e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.22 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.169106 seconds (Warm-up)
## Chain 2: 0.15076 seconds (Sampling)
## Chain 2: 0.319866 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.180245 seconds (Warm-up)
## Chain 3: 0.152475 seconds (Sampling)
## Chain 3: 0.33272 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 2e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.2 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.168738 seconds (Warm-up)
## Chain 4: 0.151746 seconds (Sampling)
## Chain 4: 0.320484 seconds (Total)
## Chain 4:summary(m3.1)## Family: binomial
## Links: mu = logit
## Formula: mc | trials(size) ~ Treatment * Time + (1 | tankID)
## Data: a (Number of observations: 23)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~tankID (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.61 0.23 0.25 1.14 1.00 1491 2108
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 3.02 0.38 2.32 3.82 1.00 1869
## Treatment2800 -1.19 0.50 -2.20 -0.22 1.00 1811
## Time80 -0.10 0.41 -0.90 0.68 1.00 2317
## Treatment2800:Time80 1.03 0.49 0.07 2.05 1.00 2332
## Tail_ESS
## Intercept 2430
## Treatment2800 2238
## Time80 2491
## Treatment2800:Time80 2226
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(marginal_effects(m3.1), ask = FALSE)## Setting the number of trials to 1 by default if not specified otherwise.
brm(mc | trials(size)~ Treatment*Time + (1|tankID), data=b,family = binomial())## Compiling the C++ model## recompiling to avoid crashing R session## Start sampling##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 3.3e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.33 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 3.51891 seconds (Warm-up)
## Chain 1: 7.55509 seconds (Sampling)
## Chain 1: 11.074 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1.4e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 4.62815 seconds (Warm-up)
## Chain 2: 8.041 seconds (Sampling)
## Chain 2: 12.6691 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.3e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 4.74891 seconds (Warm-up)
## Chain 3: 6.5455 seconds (Sampling)
## Chain 3: 11.2944 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.3e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 2.90427 seconds (Warm-up)
## Chain 4: 7.67252 seconds (Sampling)
## Chain 4: 10.5768 seconds (Total)
## Chain 4:## Warning: There were 2337 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess## Family: binomial
## Links: mu = logit
## Formula: mc | trials(size) ~ Treatment * Time + (1 | tankID)
## Data: b (Number of observations: 23)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~tankID (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.78 0.35 0.31 1.61 1.00 1077 1865
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 41.91 39.30 7.50 134.80 1.01 603
## Treatment2800 -40.04 39.31 -132.72 -5.56 1.01 605
## Time80 -38.95 39.29 -131.44 -4.52 1.01 604
## Treatment2800:Time80 39.88 39.29 5.50 132.39 1.01 603
## Tail_ESS
## Intercept 389
## Treatment2800 389
## Time80 389
## Treatment2800:Time80 389
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).brm(mc | trials(size)~ Treatment*Time + (1|tankID), data=c,family = binomial())## Compiling the C++ model## recompiling to avoid crashing R session## Start sampling##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 3.9e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 3.92987 seconds (Warm-up)
## Chain 1: 6.43746 seconds (Sampling)
## Chain 1: 10.3673 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1.7e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 3.55317 seconds (Warm-up)
## Chain 2: 7.20759 seconds (Sampling)
## Chain 2: 10.7608 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.2e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.12 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 4.67304 seconds (Warm-up)
## Chain 3: 8.58336 seconds (Sampling)
## Chain 3: 13.2564 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.4e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 5.2424 seconds (Warm-up)
## Chain 4: 4.48887 seconds (Sampling)
## Chain 4: 9.73127 seconds (Total)
## Chain 4:## Warning: There were 1749 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problems## Family: binomial
## Links: mu = logit
## Formula: mc | trials(size) ~ Treatment * Time + (1 | tankID)
## Data: c (Number of observations: 23)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Group-Level Effects:
## ~tankID (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.80 0.35 0.30 1.65 1.01 988 1880
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -41.44 33.18 -132.55 -5.36 1.01 825
## Treatment2800 43.32 33.18 7.08 133.89 1.01 828
## Time80 44.42 33.19 8.22 135.10 1.01 826
## Treatment2800:Time80 -43.49 33.20 -134.17 -7.24 1.01 828
## Tail_ESS
## Intercept 718
## Treatment2800 718
## Time80 707
## Treatment2800:Time80 732
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).OK let’s try the same thing, but with the data with 0’s for unmethylated
boxplot(mc/(mc+uc) ~ Treatment*Time, data=b, las=2)
# No random effect of tank on real data
m3 <- brm(mc | trials(size) ~ Treatment*Time, data=b, family = binomial())## Compiling the C++ model## recompiling to avoid crashing R session## Start sampling##
## SAMPLING FOR MODEL '2f81065a9d1926dcf096558e0a7da943' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 2.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.858102 seconds (Warm-up)
## Chain 1: 0.864233 seconds (Sampling)
## Chain 1: 1.72234 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '2f81065a9d1926dcf096558e0a7da943' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.671451 seconds (Warm-up)
## Chain 2: 0.790935 seconds (Sampling)
## Chain 2: 1.46239 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '2f81065a9d1926dcf096558e0a7da943' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 9e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.685373 seconds (Warm-up)
## Chain 3: 0.70976 seconds (Sampling)
## Chain 3: 1.39513 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '2f81065a9d1926dcf096558e0a7da943' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 9e-06 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.775846 seconds (Warm-up)
## Chain 4: 0.907488 seconds (Sampling)
## Chain 4: 1.68333 seconds (Total)
## Chain 4:## Warning: There were 19 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess # takes a few seconds
summary(m3)## Family: binomial
## Links: mu = logit
## Formula: mc | trials(size) ~ Treatment * Time
## Data: b (Number of observations: 23)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 39.47 33.31 7.07 133.25 1.03 113
## Treatment2800 -37.77 33.31 -131.82 -5.30 1.03 113
## Time80 -36.56 33.31 -130.57 -4.01 1.03 113
## Treatment2800:Time80 37.47 33.32 4.83 131.54 1.03 113
## Tail_ESS
## Intercept 180
## Treatment2800 180
## Time80 180
## Treatment2800:Time80 180
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m3)
# gives a similar result ot the GLM
plot(marginal_effects(m3))## Setting the number of trials to 1 by default if not specified otherwise.
# Add the random effect of tank
m3.1 <- brm(mc | trials(size)~ Treatment*Time + (1|tankID), data=b,family = binomial(), iter=5000, control = list(max_treedepth = 15))## Compiling the C++ model## recompiling to avoid crashing R session## Start sampling##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 3.3e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.33 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 1: Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 1: Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 1: Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 1: Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 1: Iteration: 2500 / 5000 [ 50%] (Warmup)
## Chain 1: Iteration: 2501 / 5000 [ 50%] (Sampling)
## Chain 1: Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 1: Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 1: Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 1: Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 1: Iteration: 5000 / 5000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 12.6971 seconds (Warm-up)
## Chain 1: 21.2541 seconds (Sampling)
## Chain 1: 33.9512 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1.4e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 2: Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 2: Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 2: Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 2: Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 2: Iteration: 2500 / 5000 [ 50%] (Warmup)
## Chain 2: Iteration: 2501 / 5000 [ 50%] (Sampling)
## Chain 2: Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 2: Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 2: Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 2: Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 2: Iteration: 5000 / 5000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 20.4341 seconds (Warm-up)
## Chain 2: 21.5739 seconds (Sampling)
## Chain 2: 42.008 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.4e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 3: Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 3: Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 3: Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 3: Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 3: Iteration: 2500 / 5000 [ 50%] (Warmup)
## Chain 3: Iteration: 2501 / 5000 [ 50%] (Sampling)
## Chain 3: Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 3: Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 3: Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 3: Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 3: Iteration: 5000 / 5000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 10.0992 seconds (Warm-up)
## Chain 3: 18.3667 seconds (Sampling)
## Chain 3: 28.4659 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.4e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 4: Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 4: Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 4: Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 4: Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 4: Iteration: 2500 / 5000 [ 50%] (Warmup)
## Chain 4: Iteration: 2501 / 5000 [ 50%] (Sampling)
## Chain 4: Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 4: Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 4: Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 4: Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 4: Iteration: 5000 / 5000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 8.24728 seconds (Warm-up)
## Chain 4: 21.548 seconds (Sampling)
## Chain 4: 29.7953 seconds (Total)
## Chain 4:summary(m3.1)## Family: binomial
## Links: mu = logit
## Formula: mc | trials(size) ~ Treatment * Time + (1 | tankID)
## Data: b (Number of observations: 23)
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup samples = 10000
##
## Group-Level Effects:
## ~tankID (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.80 0.35 0.32 1.71 1.00 2373 4452
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept 51.39 82.33 7.46 184.32 1.01 880
## Treatment2800 -49.51 82.33 -182.41 -5.58 1.01 880
## Time80 -48.41 82.33 -180.98 -4.51 1.01 874
## Treatment2800:Time80 49.33 82.33 5.44 182.20 1.01 879
## Tail_ESS
## Intercept 446
## Treatment2800 444
## Time80 448
## Treatment2800:Time80 447
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m3.1)
m3.1$data$level["contrasts"]## NULLhypothesis(m3.1, "Treatment2800 = 0", class="b", alpha=0.001)## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Treatment2800) = 0 -49.51 82.33 -1307.26 -2.45 NA
## Post.Prob Star
## 1 NA *
## ---
## 'CI': 99.8%-CI for one-sided and 99.9%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 99.9%;
## for two-sided hypotheses, the value tested against lies outside the 99.9%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.quantile(posterior_samples(m3.1)$b_Treatment2800, c(0.0001, 0.9999))## 0.01% 99.99%
## -1400.545768 -2.002076One issue with the above approach is that it requires some careful thought to do the contrasts - eg. Treatment2800 is really Treatment2800_Time09. And the Interaction is relative to Intercept + Treatment2800+Time80
OK let’s try the same thing, but with the data with 0’s for methylated Note this has some low count levels for the unmethylated section, need to fix that
boxplot(mc/(mc+uc) ~ Treatment*Time, data=c, las=2)
# No random effect of tank on real data
m3 <- brm(mc | trials(size)~ Treatment*Time, data=c,family = binomial())## Compiling the C++ model## recompiling to avoid crashing R session## Start sampling##
## SAMPLING FOR MODEL '2f81065a9d1926dcf096558e0a7da943' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 2.7e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.27 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 0.728232 seconds (Warm-up)
## Chain 1: 0.943492 seconds (Sampling)
## Chain 1: 1.67172 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '2f81065a9d1926dcf096558e0a7da943' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 9e-06 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.09 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 0.733469 seconds (Warm-up)
## Chain 2: 0.98016 seconds (Sampling)
## Chain 2: 1.71363 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '2f81065a9d1926dcf096558e0a7da943' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 7e-06 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.07 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 0.69836 seconds (Warm-up)
## Chain 3: 1.00328 seconds (Sampling)
## Chain 3: 1.70164 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '2f81065a9d1926dcf096558e0a7da943' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 3.6e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.36 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
## Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
## Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
## Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
## Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
## Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 0.787343 seconds (Warm-up)
## Chain 4: 0.815368 seconds (Sampling)
## Chain 4: 1.60271 seconds (Total)
## Chain 4:## Warning: There were 21 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded## Warning: Examine the pairs() plot to diagnose sampling problems## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#bulk-ess## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## http://mc-stan.org/misc/warnings.html#tail-ess # takes a few seconds
summary(m3)## Family: binomial
## Links: mu = logit
## Formula: mc | trials(size) ~ Treatment * Time
## Data: c (Number of observations: 23)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -45.54 37.49 -152.68 -4.63 1.01 195
## Treatment2800 47.24 37.49 6.34 154.32 1.01 195
## Time80 48.45 37.48 7.64 155.83 1.01 195
## Treatment2800:Time80 -47.54 37.48 -154.99 -6.70 1.01 195
## Tail_ESS
## Intercept 228
## Treatment2800 228
## Time80 227
## Treatment2800:Time80 228
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).pairs(m3)
# gives a similar result ot the GLM
str(marginal_effects(m3))## Setting the number of trials to 1 by default if not specified otherwise.## List of 3
## $ Treatment :'data.frame': 2 obs. of 10 variables:
## ..$ Treatment : Factor w/ 2 levels "400","2800": 1 2
## ..$ mc : num [1:2] 35.3 35.3
## ..$ size : int [1:2] 1 1
## ..$ Time : Factor w/ 2 levels "09","80": 1 1
## ..$ cond__ : Factor w/ 1 level "1": 1 1
## ..$ effect1__ : Factor w/ 2 levels "400","2800": 1 2
## ..$ estimate__: num [1:2] 5.62e-16 8.46e-01
## ..$ se__ : num [1:2] 8.34e-16 2.06e-02
## ..$ lower__ : num [1:2] 4.94e-67 8.02e-01
## ..$ upper__ : num [1:2] 0.00963 0.88318
## ..- attr(*, "effects")= chr "Treatment"
## ..- attr(*, "response")= chr "mc | trials(size)"
## ..- attr(*, "surface")= logi FALSE
## ..- attr(*, "categorical")= logi FALSE
## ..- attr(*, "ordinal")= logi FALSE
## ..- attr(*, "points")='data.frame': 12 obs. of 4 variables:
## .. ..$ Treatment: Factor w/ 2 levels "400","2800": 1 1 1 2 2 2 2 2 1 1 ...
## .. ..$ resp__ : num [1:12] 0 0 0 54 64 36 34 47 0 0 ...
## .. ..$ cond__ : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
## .. ..$ effect1__: Factor w/ 2 levels "400","2800": 1 1 1 2 2 2 2 2 1 1 ...
## $ Time :'data.frame': 2 obs. of 10 variables:
## ..$ Time : Factor w/ 2 levels "09","80": 1 2
## ..$ mc : num [1:2] 35.3 35.3
## ..$ size : int [1:2] 1 1
## ..$ Treatment : Factor w/ 2 levels "400","2800": 1 1
## ..$ cond__ : Factor w/ 1 level "1": 1 1
## ..$ effect1__ : Factor w/ 2 levels "09","80": 1 2
## ..$ estimate__: num [1:2] 5.62e-16 9.48e-01
## ..$ se__ : num [1:2] 8.34e-16 1.41e-02
## ..$ lower__ : num [1:2] 4.94e-67 9.14e-01
## ..$ upper__ : num [1:2] 0.00963 0.97155
## ..- attr(*, "effects")= chr "Time"
## ..- attr(*, "response")= chr "mc | trials(size)"
## ..- attr(*, "surface")= logi FALSE
## ..- attr(*, "categorical")= logi FALSE
## ..- attr(*, "ordinal")= logi FALSE
## ..- attr(*, "points")='data.frame': 11 obs. of 4 variables:
## .. ..$ Time : Factor w/ 2 levels "09","80": 2 2 1 1 1 2 1 1 1 2 ...
## .. ..$ resp__ : num [1:11] 36 56 0 0 0 12 0 0 0 65 ...
## .. ..$ cond__ : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
## .. ..$ effect1__: Factor w/ 2 levels "09","80": 2 2 1 1 1 2 1 1 1 2 ...
## $ Treatment:Time:'data.frame': 4 obs. of 11 variables:
## ..$ Treatment : Factor w/ 2 levels "400","2800": 1 1 2 2
## ..$ Time : Factor w/ 2 levels "09","80": 1 2 1 2
## ..$ mc : num [1:4] 35.3 35.3 35.3 35.3
## ..$ size : int [1:4] 1 1 1 1
## ..$ cond__ : Factor w/ 1 level "1": 1 1 1 1
## ..$ effect1__ : Factor w/ 2 levels "400","2800": 1 1 2 2
## ..$ effect2__ : Factor w/ 2 levels "09","80": 1 2 1 2
## ..$ estimate__: num [1:4] 5.62e-16 9.48e-01 8.46e-01 9.31e-01
## ..$ se__ : num [1:4] 8.34e-16 1.41e-02 2.06e-02 1.32e-02
## ..$ lower__ : num [1:4] 4.94e-67 9.14e-01 8.02e-01 9.00e-01
## ..$ upper__ : num [1:4] 0.00963 0.97155 0.88318 0.95451
## ..- attr(*, "effects")= chr [1:2] "Treatment" "Time"
## ..- attr(*, "response")= chr "mc | trials(size)"
## ..- attr(*, "surface")= logi FALSE
## ..- attr(*, "categorical")= logi FALSE
## ..- attr(*, "ordinal")= logi FALSE
## ..- attr(*, "points")='data.frame': 23 obs. of 6 variables:
## .. ..$ Treatment: Factor w/ 2 levels "400","2800": 1 1 1 2 1 1 2 1 2 2 ...
## .. ..$ Time : Factor w/ 2 levels "09","80": 2 2 1 2 1 1 1 2 1 2 ...
## .. ..$ resp__ : num [1:23] 36 56 0 71 0 0 54 12 64 54 ...
## .. ..$ cond__ : Factor w/ 1 level "1": 1 1 1 1 1 1 1 1 1 1 ...
## .. ..$ effect1__: Factor w/ 2 levels "400","2800": 1 1 1 2 1 1 2 1 2 2 ...
## .. ..$ effect2__: Factor w/ 2 levels "09","80": 2 2 1 2 1 1 1 2 1 2 ...
## - attr(*, "class")= chr "brmsMarginalEffects"plot(marginal_effects(m3))## Setting the number of trials to 1 by default if not specified otherwise.
# get warnings for this data
# Add the random effect of tank
m3.1 <- brm(mc | trials(size)~ Treatment*Time + (1|tankID), data=c,
family = binomial(), iter=5000, control = list(max_treedepth = 15))## Compiling the C++ model## recompiling to avoid crashing R session## Start sampling##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 3.3e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.33 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 1: Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 1: Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 1: Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 1: Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 1: Iteration: 2500 / 5000 [ 50%] (Warmup)
## Chain 1: Iteration: 2501 / 5000 [ 50%] (Sampling)
## Chain 1: Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 1: Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 1: Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 1: Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 1: Iteration: 5000 / 5000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 8.75641 seconds (Warm-up)
## Chain 1: 16.5713 seconds (Sampling)
## Chain 1: 25.3277 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1.4e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 2: Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 2: Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 2: Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 2: Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 2: Iteration: 2500 / 5000 [ 50%] (Warmup)
## Chain 2: Iteration: 2501 / 5000 [ 50%] (Sampling)
## Chain 2: Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 2: Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 2: Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 2: Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 2: Iteration: 5000 / 5000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 14.2191 seconds (Warm-up)
## Chain 2: 18.7021 seconds (Sampling)
## Chain 2: 32.9212 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.3e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 3: Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 3: Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 3: Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 3: Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 3: Iteration: 2500 / 5000 [ 50%] (Warmup)
## Chain 3: Iteration: 2501 / 5000 [ 50%] (Sampling)
## Chain 3: Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 3: Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 3: Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 3: Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 3: Iteration: 5000 / 5000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 14.0376 seconds (Warm-up)
## Chain 3: 20.2057 seconds (Sampling)
## Chain 3: 34.2434 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.3e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 4: Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 4: Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 4: Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 4: Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 4: Iteration: 2500 / 5000 [ 50%] (Warmup)
## Chain 4: Iteration: 2501 / 5000 [ 50%] (Sampling)
## Chain 4: Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 4: Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 4: Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 4: Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 4: Iteration: 5000 / 5000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 13.3362 seconds (Warm-up)
## Chain 4: 19.8597 seconds (Sampling)
## Chain 4: 33.196 seconds (Total)
## Chain 4:summary(m3.1)## Family: binomial
## Links: mu = logit
## Formula: mc | trials(size) ~ Treatment * Time + (1 | tankID)
## Data: c (Number of observations: 23)
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup samples = 10000
##
## Group-Level Effects:
## ~tankID (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.81 0.37 0.32 1.75 1.00 3071 4465
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept -43.08 37.42 -145.53 -5.14 1.00 1780
## Treatment2800 44.95 37.43 6.95 147.64 1.00 1776
## Time80 46.06 37.42 8.12 148.73 1.00 1773
## Treatment2800:Time80 -45.14 37.42 -147.82 -7.16 1.00 1772
## Tail_ESS
## Intercept 1386
## Treatment2800 1386
## Time80 1380
## Treatment2800:Time80 1387
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m3.1)
plot(marginal_effects(m3.1))## Setting the number of trials to 1 by default if not specified otherwise.
m3.1$data$Treatment["contrasts"]## [1] <NA>
## attr(,"contrasts")
## 2800
## 400 0
## 2800 1
## Levels: 400 2800m3.1$data$Time["contrasts"]## [1] <NA>
## attr(,"contrasts")
## 80
## 09 0
## 80 1
## Levels: 09 80#attr(m3.1$data, which = "terms")["term.labels"]
# (warp_em <- emmeans (m3.1, ~ Treatment))
# (warp_em <- emmeans (m3.1, ~ Time))
# (warp_em <- emmeans (m3.1, ~ Treatment | Time))
#
# cont <- contrast(warp_em, "tukey")
# cont
# (cont_posterior <- gather_emmeans_draws(cont))
# Treatment 400 and Time 09 are the reference levels (Intercept)
hypothesis(m3.1, c(
"Time80 = Treatment2800*2 + Time80 + Treatment2800:Time80", # Main effect of 400 to 2800
"Treatment2800 = Time80*2 + Treatment2800 + Treatment2800:Time80", # Main effect of day 9 to day 80
"Treatment2800:Time80=0", # Is the interaction significant
"Time80 = 0", # Is (Treatment 400 at Time 80) = (Treatment 400 and Time 09)
"Treatment2800 = 0", # Is (Treatment 2800 at Time 09) = (Treatment 400 and Time 09)
"Treatment2800 + Time80 + Treatment2800:Time80 = 0", # (Treatment 2800 at Time 80) = (Treatment 400 and Time 09)"
"Time80 = Treatment2800", # (Treatment 400 at Time 80) = (Treatment 2800 at Time 09)
"Time80 = Treatment2800 + Time80 + Treatment2800:Time80" # (Treatment 400 at Time 80) = (Treatment 2800 at Time 80)
#
), class="b", alpha=0.001)## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (Time80)-(Treatme... = 0 -44.77 37.46 -325.63 -2.81 NA
## 2 (Treatment2800)-(... = 0 -46.99 37.42 -326.70 -5.97 NA
## 3 (Treatment2800:Ti... = 0 -45.14 37.42 -325.23 -4.33 NA
## 4 (Time80) = 0 46.06 37.42 5.10 325.97 NA
## 5 (Treatment2800) = 0 44.95 37.43 3.86 325.43 NA
## 6 (Treatment2800+Ti... = 0 45.88 37.43 4.82 326.17 NA
## 7 (Time80)-(Treatme... = 0 1.11 0.65 -1.72 3.70 NA
## 8 (Time80)-(Treatme... = 0 0.18 0.67 -2.72 2.95 NA
## Post.Prob Star
## 1 NA *
## 2 NA *
## 3 NA *
## 4 NA *
## 5 NA *
## 6 NA *
## 7 NA
## 8 NA
## ---
## 'CI': 99.8%-CI for one-sided and 99.9%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 99.9%;
## for two-sided hypotheses, the value tested against lies outside the 99.9%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.plot(marginal_effects(m3.1))## Setting the number of trials to 1 by default if not specified otherwise.
# get warnings for this data for iter = 2000 and default treedepth
# get no warnings when increasing iter and treedepth
# https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
pairs(m3.1)
# in this plot want to look for iterations off the 1:1 line
# colored in red or yellow
# https://mc-stan.org/misc/warnings.htmlLet’s try the model with recoding according to levels, so we can do specific hypothesis tests a little more straightforwardly:
Results for fake data:
BASIS FOR REAL DATA HERE:
levels(c$level)## [1] "400_09" "400_80" "2800_09" "2800_80"par(mfrow=c(1,1))
boxplot(mc/size~level, data=c)
get_prior(mc | trials(size)~ level + (1|tankID), data=c,
family = binomial())## prior class coef group resp dpar nlpar bound
## 1 b
## 2 b level2800_09
## 3 b level2800_80
## 4 b level400_80
## 5 student_t(3, 0, 10) Intercept
## 6 student_t(3, 0, 10) sd
## 7 sd tankID
## 8 sd Intercept tankIDm4.1 <- brm(mc | trials(size)~ level + (1|tankID), data=c,
family = binomial(), iter=5000, control = list(max_treedepth = 15, adapt_delta=0.99))## Compiling the C++ model## recompiling to avoid crashing R session## Start sampling##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 1).
## Chain 1:
## Chain 1: Gradient evaluation took 3.6e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.36 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1:
## Chain 1:
## Chain 1: Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 1: Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 1: Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 1: Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 1: Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 1: Iteration: 2500 / 5000 [ 50%] (Warmup)
## Chain 1: Iteration: 2501 / 5000 [ 50%] (Sampling)
## Chain 1: Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 1: Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 1: Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 1: Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 1: Iteration: 5000 / 5000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 39.4546 seconds (Warm-up)
## Chain 1: 31.5092 seconds (Sampling)
## Chain 1: 70.9638 seconds (Total)
## Chain 1:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 2).
## Chain 2:
## Chain 2: Gradient evaluation took 1.4e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2:
## Chain 2:
## Chain 2: Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 2: Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 2: Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 2: Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 2: Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 2: Iteration: 2500 / 5000 [ 50%] (Warmup)
## Chain 2: Iteration: 2501 / 5000 [ 50%] (Sampling)
## Chain 2: Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 2: Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 2: Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 2: Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 2: Iteration: 5000 / 5000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 34.8664 seconds (Warm-up)
## Chain 2: 35.7707 seconds (Sampling)
## Chain 2: 70.6371 seconds (Total)
## Chain 2:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 3).
## Chain 3:
## Chain 3: Gradient evaluation took 1.4e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3:
## Chain 3:
## Chain 3: Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 3: Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 3: Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 3: Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 3: Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 3: Iteration: 2500 / 5000 [ 50%] (Warmup)
## Chain 3: Iteration: 2501 / 5000 [ 50%] (Sampling)
## Chain 3: Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 3: Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 3: Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 3: Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 3: Iteration: 5000 / 5000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 29.0224 seconds (Warm-up)
## Chain 3: 60.5825 seconds (Sampling)
## Chain 3: 89.6049 seconds (Total)
## Chain 3:
##
## SAMPLING FOR MODEL '5ebef1e8424272f9f2d56ea1a767a123' NOW (CHAIN 4).
## Chain 4:
## Chain 4: Gradient evaluation took 1.3e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4:
## Chain 4:
## Chain 4: Iteration: 1 / 5000 [ 0%] (Warmup)
## Chain 4: Iteration: 500 / 5000 [ 10%] (Warmup)
## Chain 4: Iteration: 1000 / 5000 [ 20%] (Warmup)
## Chain 4: Iteration: 1500 / 5000 [ 30%] (Warmup)
## Chain 4: Iteration: 2000 / 5000 [ 40%] (Warmup)
## Chain 4: Iteration: 2500 / 5000 [ 50%] (Warmup)
## Chain 4: Iteration: 2501 / 5000 [ 50%] (Sampling)
## Chain 4: Iteration: 3000 / 5000 [ 60%] (Sampling)
## Chain 4: Iteration: 3500 / 5000 [ 70%] (Sampling)
## Chain 4: Iteration: 4000 / 5000 [ 80%] (Sampling)
## Chain 4: Iteration: 4500 / 5000 [ 90%] (Sampling)
## Chain 4: Iteration: 5000 / 5000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 25.5027 seconds (Warm-up)
## Chain 4: 30.739 seconds (Sampling)
## Chain 4: 56.2417 seconds (Total)
## Chain 4:summary(m4.1)## Family: binomial
## Links: mu = logit
## Formula: mc | trials(size) ~ level + (1 | tankID)
## Data: c (Number of observations: 23)
## Samples: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
## total post-warmup samples = 10000
##
## Group-Level Effects:
## ~tankID (Number of levels: 12)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 0.80 0.37 0.31 1.70 1.00 2384 3968
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept -44.70 41.72 -159.76 -5.20 1.00 1419 723
## level400_80 47.69 41.73 8.19 162.74 1.00 1420 728
## level2800_09 46.57 41.73 6.99 161.68 1.00 1423 728
## level2800_80 47.50 41.73 7.97 162.43 1.00 1423 727
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).plot(m4.1)
m4.1$data$level["contrasts"]## [1] <NA>
## attr(,"contrasts")
## 400_80 2800_09 2800_80
## 400_09 0 0 0
## 400_80 1 0 0
## 2800_09 0 1 0
## 2800_80 0 0 1
## Levels: 400_09 400_80 2800_09 2800_80str(posterior_samples(m4.1))## 'data.frame': 10000 obs. of 18 variables:
## $ b_Intercept : num -52.3 -42.9 -37.6 -11.9 -53 ...
## $ b_level400_80 : num 55.1 45.4 41.4 15.3 55.2 ...
## $ b_level2800_09 : num 53.9 44.8 39.1 13.5 55.2 ...
## $ b_level2800_80 : num 55 45.4 40.6 15.3 56.1 ...
## $ sd_tankID__Intercept : num 0.0842 0.0795 0.0635 0.3977 0.8953 ...
## $ r_tankID[1,Intercept] : num -0.0937 0.0855 0.0333 0.2038 2.0906 ...
## $ r_tankID[2,Intercept] : num 0.1025 -0.0492 -0.0485 0.1975 -0.3963 ...
## $ r_tankID[3,Intercept] : num -0.04354 0.05142 0.00843 0.63892 -2.0433 ...
## $ r_tankID[4,Intercept] : num -0.0592 0.0846 0.035 0.0136 0.6248 ...
## $ r_tankID[5,Intercept] : num -0.0629 0.0168 0.0141 -0.12 0.1004 ...
## $ r_tankID[6,Intercept] : num 0.04353 -0.12189 -0.00554 -0.3019 -1.20841 ...
## $ r_tankID[7,Intercept] : num 0.0434 -0.0919 -0.0698 0.5376 0.5407 ...
## $ r_tankID[8,Intercept] : num 0.0311 -0.04 0.0217 -0.618 -1.339 ...
## $ r_tankID[9,Intercept] : num 0.1288 -0.1769 -0.0883 1.0328 0.7092 ...
## $ r_tankID[10,Intercept]: num 0.0759 -0.0224 -0.0135 0.7691 0.0797 ...
## $ r_tankID[11,Intercept]: num -0.1116 0.0438 0.0273 -0.6215 -1.3675 ...
## $ r_tankID[12,Intercept]: num -0.0316 0.0233 -0.0464 0.2665 0.3883 ...
## $ lp__ : num -77.1 -82.3 -84.4 -77.5 -74.9 ...p<-plot(marginal_effects(m4.1))## Setting the number of trials to 1 by default if not specified otherwise.
round(p$level$data$estimate,3)## [1] 0.000 0.951 0.865 0.942round(p$level$data$upper__,3)## [1] 0.006 0.983 0.938 0.976round(p$level$data$lower__,3)## [1] 0.000 0.885 0.753 0.883marginal_effects(m4.1)## Setting the number of trials to 1 by default if not specified otherwise.
alph = 0.0001/2
hypothesis(m4.1, c(
"level2800_09 + 2*Intercept= 2*Intercept + level2800_80 + level400_80",# Is Day 9 is different from Day 80
"level400_80 + 2*Intercept= 2*Intercept + level2800_80 + level2800_09", # Is 400 different from 2800
"Intercept-(Intercept + level400_80) = (Intercept + level2800_09) - (Intercept + level2800_80)", # Is interaction significant
"Intercept-(Intercept + level2800_09) = (Intercept + level400_80) - (Intercept + level2800_80)", # Is interaction significant (should be same as last row)
"level2800_80>0", # Is 2800_80 > 400_09
"level2800_80<0", # Is 2800_80 > 400_09
"level2800_09<0",
"level2800_09>0",
"level400_80>0",
"level400_80<0",
"level2800_80=level2800_09",
"level2800_80=level400_80",
"level2800_09=level400_80"
), class="b", alpha=alph)## Hypothesis Tests for class b:
## Hypothesis Estimate Est.Error CI.Lower CI.Upper
## 1 (level2800_09+2*I... = 0 -48.62 41.73 -491.92 -5.18
## 2 (level400_80+2*In... = 0 -46.38 41.73 -491.41 -1.64
## 3 (Intercept-(Inter... = 0 -46.76 41.73 -489.85 -3.42
## 4 (Intercept-(Inter... = 0 -46.76 41.73 -489.85 -3.42
## 5 (level2800_80) > 0 47.50 41.73 4.57 485.78
## 6 (level2800_80) < 0 47.50 41.73 4.57 485.78
## 7 (level2800_09) < 0 46.57 41.73 3.73 484.89
## 8 (level2800_09) > 0 46.57 41.73 3.73 484.89
## 9 (level400_80) > 0 47.69 41.73 4.47 485.37
## 10 (level400_80) < 0 47.69 41.73 4.47 485.37
## 11 (level2800_80)-(l... = 0 0.93 0.27 0.00 2.32
## 12 (level2800_80)-(l... = 0 -0.19 0.66 -5.23 2.96
## 13 (level2800_09)-(l... = 0 -1.12 0.64 -6.07 1.98
## Evid.Ratio Post.Prob Star
## 1 NA NA *
## 2 NA NA *
## 3 NA NA *
## 4 NA NA *
## 5 Inf 1 *
## 6 0 0
## 7 0 0
## 8 Inf 1 *
## 9 Inf 1 *
## 10 0 0
## 11 NA NA
## 12 NA NA
## 13 NA NA
## ---
## 'CI': 99.99%-CI for one-sided and 99.995%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 99.995%;
## for two-sided hypotheses, the value tested against lies outside the 99.995%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.Results for real data:
Note that I had this error message: Warning messages: 1: There were 1 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help.
And as described in the paper I increased adapt_delta to 0.99
Plan: Use last block of code as a baseline:
Read up more on how to set priors and think about reasonalbe priors to set. The priors for the binomial are set with the student_t distribution.
DONE. Order the levels so that 400_09 is the first level and recode the hypothesis tests accordingly. (Check that my coding of the hypothesis tests are correct first! - Need to understand “contrasts” better). Then go over re-coding together to make sure it is correct. Understand the contrasts to make sure hypothesis test is correct.
DONE. Read up more on the hypothesis tests and make sure we understand it.
?hypothesisis a good place to start, and read Kass and Rafferty https://www.tandfonline.com/doi/abs/10.1080/01621459.1995.10476572 The function help says to avoid using default priors. Note that evidence ratio can be calculated for one-sided hypotheses but not two-sided hypotheses, and the evidence ratio can be interpreted as strength of evidence of an effect. I think we want to do one-sided hypotheses with alpha/2 so we can store the Evidence Ratio, which according to Kass and Raftery log10(Evidence Ratio) > 1 is “strong evidence” and log10(Evidence Ratio) > 2 is “decisive evidence”Decide on alpha. This is Bayesian, so no way to correct for multiple tests. Here the help on
hypothesisis also very good. There are 10,000 posterior samples, soalpha=0.0001means that only 1 of the posterior samples will lie outside the credible interval.Write up the methods for the paper. Here is a start:
Methylation data is binomially distributed (e.g., read counts for methylated or unmethylated) and in our case, we wanted to make sure to account for the random effect of tank in the experiment (see results). (Current packages such as methylkit and MACAU do not have the ability to include random effects). For some loci, we found with preliminary analyses that some approaches to generalized linear modeling failed to converge (when parameters were estimated via restricted maximum likelihood as implemented in lme4, or when the posterior distributions of the parameters were estimated via MCMC with Gibbs sampling as implemented in MCMCglmm). We traced this issue of convergence to loci in which all individuals within a treatment combination had 100% or 0% methylation, and with diverse levels of methylation in other treatments, making them biologically interesting and desireable to model correctly. As pointed out by Burker, the approaches used by MCMCglmm can be slow to converge in high-dimensional models with correlated parameters and may depend on conjugate priors.
We therefore took a robust approach to analyzing metylation data with Bayesian Regression Models (BRMs) in which the sampling on the posterior distributions on the parameters was accomplished with the No-U-Turn-Sampler (NUTS). This algorithm converges much more quickly for high-dimensional models regardless of whether the priors are conjugate or not. To avoid divergent transitions, to ensure that transitions after warmup did not exceed the maximum treedepth, and to ensure the effective sample size fo the bulk and tail was sufficient, we increased adapt_delta to 0.99, max_treedepth to 15, the number of transitions to 5000, respectively. We modeled percent methylation as a binomial response variable with treatment level (treatment and timepoint) as a fixed/population effect and tank as a random/group effect. We then used Bayes Factors to test three hypotheses: (i) main effect of treatment (400 ppm = 2800 ppm), (ii) main effect of time (day 9 = day 80), and (iii) an interaction between time and treatment. Bayes Factors were calculated with the hypothesis function from brms. Note that in the Bayesian framework, P-values are not calculated and therefore there is no analog to a false discovery rate correction. Hypothesis tests with log10(Bayes Factors) > 2 and non-overlapping posterior distributions (corresponding to 10,000 posterior samples) were taken to be strong evidence of an effect. For loci that had a significant interaction between time and treatment, we conduced an additional set of post-hoc hypothesis tests of individual treatment levels using the same criteria as above.
With n posterior draws, the number of false positives would be expected to be 1/n*(number of loci)*(3 tests per loci) + 1/n*(number of loci with significant interaction)*(6 comparisons).
ALAN. Write a function to test for main effects and an interaction, if the interaction is significant than also test for overlap of CI for pairwise params. For signifcant params store relevant outputs. Check that all output is there. Test on top 10 loci and 10 random loci to make sure everything works
Save the “Estimate” and the “Star” (h\(hypothesis\)Star=="*“) from each hypothesis test, but be sure to note”alpha" in the results (e.g.
h_2800.80E2800.09_a0.0001)ALAN. Add logical to the function to output the plot if significant, but make sure to output a consistent plot for each locus: https://github.com/paul-buerkner/brms/issues/59
Pseudoparallize code to run (n loci)/(68 cores) on each core
Other notes:
Note that I do not think we have a zero-inflated model. This model assumes that the sample is a “mixture” of two sorts of individuals: one group whose counts are generated by the standard binomial regression model, and another group (call them the absolute zero group) who have zero probability of a count greater than 0.
In our case, we may have 0% or 100% methylation, but the counts are greater than 0 in at least one of the categories.